# Gauss积分和自适应Gauss积分的简要实现
# 真正运用时，可以选择现成的计算库

import sympy
import sympy.integrals.quadrature
import numpy as np

def gauss_int(f,a,b,npoint = 3):
    a = float(a)
    b = float(b)
    npoint = int(npoint)

    if a >= b or npoint <= 1:
        raise ValueError('wrong input')

    xi, wi = sympy.integrals.quadrature.gauss_legendre(npoint,5)
    xi, wi = np.array(xi,dtype=np.float64)[:,None],np.array(wi,dtype=np.float64)[:,None]

    I = (b-a)/2*wi.T@f((a+b)/2+(b-a)/2*xi);
    return I.item()

def adaptive_gauss_int(f,a,b,npoint = 3,err = 0.001,depth = 0):
    a = float(a)
    b = float(b)
    err = float(err)
    npoint = int(npoint)
    depth = int(depth)

    c = (a+b)/2
    I =  gauss_int(f,a,b,npoint)
    I1 = gauss_int(f,a,c,npoint)
    I2 = gauss_int(f,c,b,npoint)

    if depth > 5:
        print('not converged')
        return I1 + I2
    
    if np.abs(I-I1-I2) > err:
        new_I = adaptive_gauss_int(f,a,c,npoint,err,depth+1) + adaptive_gauss_int(f,c,b,npoint,err,depth+1)
        return new_I
    else:
        return I1 + I2

x = sympy.symbols('x')
expr = sympy.cos(x)
func = sympy.lambdify((x),expr,'numpy')

I = sympy.integrate(expr,(x,0,20)).evalf()
print(I)

I = gauss_int(func,0,20)
print(I)

I = adaptive_gauss_int(func,0,20)
print(I)
